#########################################################################################
####################################  FIGURE 6   ########################################
#######  Simulated voting rates, by condition and intensity of discrimination  ##########
#########################################################################################

library(pacman)
p_load(tidyverse, ggplot2, msm, sandwich, here, xtable)
set.seed(12345)
rm(list = ls())


#########################################################################################

load(file = "simulations/simulations_output_2.Rda")

# predict discrimination rate from observed X in inverted function
getX <- function(y, model) {
  b <- coef(model)[1]
  c <- coef(model)[2]
  d <- coef(model)[3]
  x <- -(sqrt(-4*b*d+ c^2 + 4*d*y) + c)/(2*d)
  return(x)
}

outcomes <- factor(c("Team-building influence",
                     "Global assessment influence",
                     "Team-building\nincentivized spokesperson",
                     "Global assessment,\nbefore leader",
                     "Global assessment,\nmale leader",
                     "Global assessment,\nfemale leader",
                     "Team-building influence,\nbefore leader",
                     "Team-building influence,\nmale leader",
                     "Team-building influence,\nfemale leader",
                     "Global influence,\nMonth 1",
                     "Global influence,\nMonth 2",
                     "Global influence,\nMonth 3",
                     "Global influence,\nMonth 4"), levels =
                     c("Team-building influence",
                       "Global assessment influence",
                       "Team-building\nincentivized spokesperson",
                       "Global assessment,\nbefore leader",
                       "Global assessment,\nmale leader",
                       "Global assessment,\nfemale leader",
                       "Team-building influence,\nbefore leader",
                       "Team-building influence,\nmale leader",
                       "Team-building influence,\nfemale leader",
                       "Global influence,\nMonth 1",
                       "Global influence,\nMonth 2",
                       "Global influence,\nMonth 3",
                       "Global influence,\nMonth 4"))

# means -- taken from graphs of figure 1 and related
# survey measures are from "surveyminfl_bymonth_forgraphing.csv"

# Data points from Figure 1:
task_influence <- read.csv("figure_1/intermediate_data/combined_groupgenderminfl_ratio_ravg_forgraphing.csv")
global_influence <- read.csv("figure_1/intermediate_data/combined_groupfemvotesminfl_avgpre_ratio_forgraphing.csv")
spokesperson <- read.csv("figure_1/intermediate_data/combined_genderspokes_ratio_ravg_forgraphing.csv")

# Data points from Figure 3:
minfl_ratio_r1 <- read.csv("figure_1/intermediate_data/study2_groupgenderminfl_ratio_r1_forgraphing.csv")
minfl_ratio_r2_fem <- read.csv("figure_3/intermediate_data/study2_newfig_groupgenderminfl_ratio_r2_fem_forgraphing.csv")
minfl_ratio_r2_mal <- read.csv("figure_3/intermediate_data/study2_newfig_groupgenderminfl_ratio_r2_male_forgraphing.csv")

# Data points from Figure 4:
surveyminfl_bymonth <- read.csv("figure_4/intermediate_data/surveyminfl_bymonth_forgraphing.csv")

# Data points from global assesment data:
femvotesminfl_pre <- read.csv("figure_6/intermediate_data/study2_newfig_groupfemvotesminfl_avgpre_ratio_forgraphing.csv")
femvotesminfl_fem <- read.csv("figure_6/intermediate_data/study2_newfig_groupfemvotesminfl_avgpost_ratio_fem_forgraphing.csv")
femvotesminfl_male <- read.csv("figure_6/intermediate_data/study2_newfig_groupfemvotesminfl_avgpost_ratio_male_forgraphing.csv")


femMajValues  <- c(task_influence$v2[1], global_influence$v2[1], spokesperson$v2[1], # From Figure 1
                   femvotesminfl_pre$v2[1], femvotesminfl_male$v2[1], femvotesminfl_fem$v2[1], # From Global Assesment data.
                   minfl_ratio_r1$v2[1], # From Figure 1
                   minfl_ratio_r2_mal$v2[1], minfl_ratio_r2_fem$v6[1], # From Figure 3
                   surveyminfl_bymonth$estimate[2], surveyminfl_bymonth$estimate[4], # From Figure 4
                   surveyminfl_bymonth$estimate[6], surveyminfl_bymonth$estimate[8]) # From Figure 4

maleMajValues <- c(task_influence$v2[2], global_influence$v2[2], spokesperson$v2[2], # From Figure 1
                   femvotesminfl_pre$v2[2], femvotesminfl_male$v2[2], femvotesminfl_fem$v2[2],  # From Global Assesment data.
                   minfl_ratio_r1$v2[2], # From Figure 1
                   minfl_ratio_r2_mal$v2[2], minfl_ratio_r2_fem$v6[2], # From Figure 3
                   surveyminfl_bymonth$estimate[1], surveyminfl_bymonth$estimate[3], 
                   surveyminfl_bymonth$estimate[5], surveyminfl_bymonth$estimate[7]) # From Figure 4



#femMajValues  <- c(.793, .856, 0.772, 0.856, 0.7435, 0.9615, .73, .641, .854, 0.824, 0.891, 0.856, 0.896)
#maleMajValues <- c(.447, .606, 0.499, 0.6965, 0.627, 0.899, .438, .419, .714, 0.376, 0.673, 0.731, 0.91)



femMajDiscrim <- getX(femMajValues, femaleMajorityTrend)
maleMajDiscrim <- getX(maleMajValues, maleMajorityTrend)

pointsFrame <- as.data.frame(cbind(femMajDiscrim, femMajValues, maleMajDiscrim, maleMajValues))
pointsFrame <- cbind(pointsFrame, outcomes)


loopOut.df.aggregated$predictionFemMaj <- getX(loopOut.df.aggregated$FvoteRatioBest3F, femaleMajorityTrend)
loopOut.df.aggregated$predictionMaleMaj <- getX(loopOut.df.aggregated$FvoteRatioBest1F, maleMajorityTrend)

# set graph option parameters
colorFemaleMaj<-"#E69F00"
colorMaleMaj<- "#56B4E9"
gphLinewidth <- 1.5
gphPointSize <- 5
gphAxisLabFontSize <- 28

# outcomes
subpointsFrame<- pointsFrame[1:3,]
ggplot(data = loopOut.df.aggregated) +
  geom_line(aes(x = predictionMaleMaj, y = FvoteRatioBest1F, color=colorMaleMaj ), linewidth = gphLinewidth) +
  geom_line(aes(x = predictionFemMaj, y = FvoteRatioBest3F, color=colorFemaleMaj), linewidth = gphLinewidth) +
  geom_point(data = subpointsFrame, aes(x = femMajDiscrim, y = femMajValues, shape = outcomes), color = colorFemaleMaj, size = gphPointSize) +
  geom_point(data = subpointsFrame, aes(x = maleMajDiscrim, y = maleMajValues, shape = outcomes), color = colorMaleMaj, size = gphPointSize) +
  ylab("votes per woman") + xlab("Performance penalty for women (SD)") +
  geom_hline(yintercept = 1, linetype = 'dashed', color="black") +
  scale_color_manual(values=c(colorMaleMaj, colorFemaleMaj),
                     labels = c("Male-majority vote series",
                                "Female-majority vote series")) +
  theme(legend.title=element_blank(), legend.position = c(.8, .85)) +
  theme(text=element_text(size=12),
        axis.title=element_text(size=gphAxisLabFontSize),
        axis.text=element_text(size=gphAxisLabFontSize-6),
        plot.title = element_text(hjust=.5),
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        panel.background = element_blank(),
        axis.title.x=element_blank(),
        axis.line = element_line(colour="black"))+
  guides(shape = guide_legend(override.aes = list(color = "grey"))) +
  scale_y_continuous(breaks=seq(0,1.0,.2), expand = c(0,0), limits = c(0, 1.2)) +   scale_x_continuous(breaks=seq(-.10,.8,.10), expand = c(0,0), limits = c(-.2,.9))


ggsave(file = 'figure_6/figure_6a.pdf', width = 12, height = 7, units = "in")


# Leadership, global influence
subpointsFrame<- pointsFrame[4:6,]
ggplot(data = loopOut.df.aggregated) +
  geom_line(aes(x = predictionMaleMaj, y = FvoteRatioBest1F, color=colorMaleMaj ), linewidth = gphLinewidth) +
  geom_line(aes(x = predictionFemMaj, y = FvoteRatioBest3F, color=colorFemaleMaj), linewidth = gphLinewidth) +
  geom_point(data = subpointsFrame, aes(x = femMajDiscrim, y = femMajValues, shape = outcomes), color = colorFemaleMaj, size = gphPointSize) +
  geom_point(data = subpointsFrame, aes(x = maleMajDiscrim, y = maleMajValues, shape = outcomes), color = colorMaleMaj, size = gphPointSize) +
  ylab("votes per woman") + xlab("Performance penalty for women (SD)") +
  geom_hline(yintercept = 1, linetype = 'dashed', color="black") +
  scale_color_manual(values=c(colorMaleMaj, colorFemaleMaj),
                     labels = c("Male-majority vote series",
                                "Female-majority vote series")) +
  theme(legend.title=element_blank(), legend.position = c(.8, .85)) +
  theme(text=element_text(size=12),
        axis.title=element_text(size=gphAxisLabFontSize),
        axis.text=element_text(size=gphAxisLabFontSize-6),
        plot.title = element_text(hjust=.5),
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        panel.background = element_blank(),
        axis.line = element_line(colour="black"))+
  guides(shape = guide_legend(override.aes = list(color = "grey"))) +
  scale_y_continuous(breaks=seq(0,1.0,.2), expand = c(0,0), limits = c(0, 1.2)) +   scale_x_continuous(breaks=seq(-.10,.8,.10), expand = c(0,0), limits = c(-.2,.9))

ggsave( file = 'figure_6/figure_6c.pdf', width = 12, height = 7, units = "in")


# leadership task
subpointsFrame<- pointsFrame[7:9,]
ggplot(data = loopOut.df.aggregated) +
  geom_line(aes(x = predictionMaleMaj, y = FvoteRatioBest1F, color=colorMaleMaj ), linewidth = gphLinewidth) +
  geom_line(aes(x = predictionFemMaj, y = FvoteRatioBest3F, color=colorFemaleMaj), linewidth = gphLinewidth) +
  geom_point(data = subpointsFrame, aes(x = femMajDiscrim, y = femMajValues, shape = outcomes), color = colorFemaleMaj, size = gphPointSize) +
  geom_point(data = subpointsFrame, aes(x = maleMajDiscrim, y = maleMajValues, shape = outcomes), color = colorMaleMaj, size = gphPointSize) +
  ylab("votes per woman") + xlab("Performance penalty for women (SD)") +
  geom_hline(yintercept = 1, linetype = 'dashed', color="black") +
  scale_color_manual(values=c(colorMaleMaj, colorFemaleMaj),
                     labels = c("Male-majority vote series",
                                "Female-majority vote series")) +
  theme(legend.title=element_blank(), legend.position = c(.8, .85)) +
  theme(text=element_text(size=12),
        axis.title=element_text(size=gphAxisLabFontSize),
        axis.text=element_text(size=gphAxisLabFontSize-10),
        plot.title = element_text(hjust=.5),
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        panel.background = element_blank(),
        axis.title.x=element_blank(),
        axis.title.y=element_blank(),
        axis.line = element_line(colour="black"))+
  guides(shape = guide_legend(override.aes = list(color = "grey"))) +
  scale_y_continuous(breaks=seq(0,1.0,.2), expand = c(0,0), limits = c(0, 1.2)) +   scale_x_continuous(breaks=seq(-.10,.8,.10), expand = c(0,0), limits = c(-.2,.9))

ggsave(file = 'figure_6/figure_6b.pdf', width = 12, height = 7, units = "in")


subpointsFrame<- pointsFrame[10:13,]
subpointsFrame$monthLabel <- seq(1,4,1)

ggplot(data = loopOut.df.aggregated) +
  geom_line(aes(x = predictionMaleMaj, y = FvoteRatioBest1F, color=colorMaleMaj ), linewidth = gphLinewidth) +
  geom_line(aes(x = predictionFemMaj, y = FvoteRatioBest3F, color=colorFemaleMaj), linewidth = gphLinewidth) +
  geom_text(data = subpointsFrame, aes(x = femMajDiscrim, y = femMajValues), label = subpointsFrame$monthLabel, size = gphPointSize+1) +
  geom_text(data = subpointsFrame, aes(x = maleMajDiscrim, y = maleMajValues), label = subpointsFrame$monthLabel, size = gphPointSize+1) +
  ylab("votes per woman") + xlab("Performance penalty for women (SD)") +
  geom_hline(yintercept = 1, linetype = 'dashed', color="black") +
  scale_color_manual(values=c(colorMaleMaj, colorFemaleMaj),
                     labels = c("Male-majority vote series",
                                "Female-majority vote series")) +
  theme(legend.title=element_blank(), legend.position = c(.8, .85)) +
  theme(text=element_text(size=12),
        axis.title=element_text(size=gphAxisLabFontSize),
        axis.text=element_text(size=gphAxisLabFontSize-6),
        plot.title = element_text(hjust=.5),
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        panel.background = element_blank(),
        axis.title.y=element_blank(),
        axis.line = element_line(colour="black"))+
  guides(shape = guide_legend(override.aes = list(color = "grey"))) +
  scale_y_continuous(breaks=seq(0,1.0,.2), expand = c(0,0), limits = c(0, 1.2)) +   scale_x_continuous(breaks=seq(-.10,.8,.10), expand = c(0,0), limits = c(-.2,.9))

ggsave(file = 'figure_6/figure_6d.pdf', width = 12, height = 7, units = "in")